C===============================================================
C  Make the Green's function on X mesh
      subroutine gf18(ln,eproj,eta,ucentr,cntfug,gridx,nx,
     >   jdouble,njdouble,regcut,expcut,reg,beta,jstart,jstop,phase)
c
      implicit real*8 (a-h, o-z)
c
      dimension ucentr(nx),cntfug(nx),gridx(nx),jdouble(njdouble),
     >   reg(nx),beta(nx)
      complex*16 phase
      complex*16 xx,eta1,zlmin,cfc(1),cgc(1),cfcp(1),
     >   cgcp(1),sig(1)

      call reglr8(ln,eproj,eta,ucentr,cntfug,gridx,nx,
     >   jdouble,njdouble,regcut,expcut,reg,jstart,jstop,phase)

      xlmin = dble(ln)
      zlmin = dcmlpx(xlmin,0d0)
      ifail = 1
      nl = 1
      
C  Now to find the iregular solution. We take care so as not to start
C  at values that are too small. Use routine COULCC to get the G function
      j = jstop
      if (eproj.lt.0d0) then
         phase = phase/(0d0,1d0)
         mode1 = 12
         kfn = 0
         eta1 = eta/(0d0,1d0)
         do while (j.ge.jstop-1)
            xx = cdsqrt(dcmplx(eproj)) * dcmplx(gridx(j))
            call coulcc(xx,eta1,zlmin,nl,cfc,cgc,cfcp,cgcp,sig,mode1,
     >         kfn,ifail)
            beta(j) = abs(cgc(1))
            j = j-1
         end do 
         phase = phase*cgc(1)/abs(cgc(1))
         rph = dreal(phase)
         aph = dble(dimag(phase))
         if (abs(aph).gt.1d-5) then
            print*,'Problems in making Green''s functions'
            print*,'Imaginary part of phase is too big, expect 0',aph
         end if 
      else if (eproj.eq.0d0) then
         do while (j.ge.jstop-1)
            if (eta.eq.0d0) then
               beta(j) = gridx(j)**(-ln)
            else
C  Here to evaluate the Coulomb G.F. at zero energy
               xx = dcmplx(dsqrt(8d0*gridx(j)))
               mode1 = 2
               kfn = 2
               xlmin = dble(2*ln+1)
               zlmin = dcmlpx(xlmin, 0d0)
               call coulcc(xx,eta1,zlmin,nl,cfc,cgc,cfcp,cgcp,sig,
     >            mode1,kfn,ifail)
               beta(j) = -dreal(cgc(1))*dsqrt(gridx(j))
            end if
            j = j-1
         end do 
      else
         jj = jstop
         do while (j.ge.jstop-1.or.abs(beta(jj)) .lt. 0.1d0)
            xx = cdsqrt(dcmplx(eproj)) * dcmplx(gridx(j))
            eta1 = dcmplx(eta,0d0)
            mode1 = 2
            kfn = 0
            call coulcc(xx,eta1,zlmin,nl,cfc,cgc,cfcp,cgcp,sig,mode1,
     >         kfn,ifail)
            beta(j) = dreal(cgc(1))
            jj = j
            j = j-1
         end do
      end if 
      
      j = j+2
      s2 = beta(j-1)
      s3 = beta(j)
      call nmrvb8(ln,eproj,ucentr,cntfug,gridx,nx,
     >   jdouble,njdouble,s2,s3,beta,jstart,j)
      t2 = 1d0/dreal(phase)
      t = t2 * dble(imag(phase))
      if (abs(t).gt.1e-4) then
         DO J = JSTART,JSTOP
            BETA(J) = beta(j)*t2-t*reg(j)
         end do
      else
C  T2 may be -1d0         
         DO J = JSTART,JSTOP
            BETA(J) = beta(j)*t2
         end do
      end if 
      RETURN
      END
